#ifndef ARRAYVECTORCLASS_H
#define ARRAYVECTORCLASS_H

#include <string>
#include <iostream>
#include <ostream>
#include <fstream>
#include <cstdlib>

using std::ostream;
using std::string;
using std::istream;
// templates for arrays of dimension 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9.
// test of indices not included.

// Dimensions: d0 , d1 , d2 , d3 , d4 , d5 , d6 , d7 , d8 for each dimension.
// The dimension of one component is returned by dimension: A.dimension (4) = d4 for example
// Constructors and allocators of N dimension arrays are given the N dimensions of the arrays as inputs.
// One has copy constructors from class array or a standard pointer array , in which case the dimension is 1.
// The unused dimensions are put to zero.
// The total number of elements is dimension_table , product of the used dimensions. It is given by A.dimension_total ()
//
// One can test if an array is allocated with is_it_filled
//
// The array of elements is a standard allocated pointer to C template type inside the class.
// An element is accessed from its N components with A(i0 , i1 , i2) for example.
// The function index_determine sends the one-dimension index i associated to the N indices: A.index_determine(i0 , i1 , i2) = i2 + d2 * (i1 + d1 * i0) for example.
// An element is accessed from its one-dimensional index with A[i] for example , with i = A.index_determine (i0 , i1 , i2) .
//
// One can allocate and copy an array of user-defined classes with allocate_fill_object_elements.
// One can deallocate the allocated objects of the array without deallocating the array itself with deallocate_object_elements.
//
// Standard operator overloadings of the form = , + , - , * , / and -= , += , *= , /= can be used.
// Comparison operators of form == , != can be used.
//
// One can copy and read an array from disk with copy_disk (output.dat) and read_disk (input.dat) for example. Files are not formatted.
//
// Sum , product , min and max of array elements can be calculated.
//
// The used memory of the class can be calculated with used_memory_calc.
// A recursive search is done if C is a class. used_memory_calc must be defined for all classes in this recursive calculation.
// used_memory_calc is defined for all basic types.
// It is in Mb.
//
// One can sort the array with quick_sort (general) , quick_sort_by_real_part (sorting by the real parts of elements) , quick_sort_by_abs (sorting by the moduli of elements)
//
// One can reverse the order of elements with reverse_order
//
// One can assign the content of an array in an array of different dimension with assign
//
// One can find the index of a value closest to a given value with closest_value_index_determine
//
// One can check if an array is equally spaced with is_it_equally_spaced
//
// Standard MPI routines are implemented (see MPI_helper)



// The next routine does nothing in the usual mode. But it would print "Test mode for arrays and vectors indices" on screen with testing mode to remind that one is in the array-test mode , which is much slower than the usual mode.

template <class C>
void print_array_vector_test_status () {}

template <class C>
class array
{
public:

  explicit array () : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL) {}

  explicit array (const unsigned int d0_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c ,
		  const unsigned int d4_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c , d4_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c ,
		  const unsigned int d4_c ,
		  const unsigned int d5_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c , d4_c , d5_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c ,
		  const unsigned int d4_c ,
		  const unsigned int d5_c ,
		  const unsigned int d6_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c , d4_c , d5_c , d6_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c ,
		  const unsigned int d4_c ,
		  const unsigned int d5_c ,
		  const unsigned int d6_c ,
		  const unsigned int d7_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c , d4_c , d5_c , d6_c , d7_c);
  }

  explicit array (
		  const unsigned int d0_c ,
		  const unsigned int d1_c ,
		  const unsigned int d2_c ,
		  const unsigned int d3_c ,
		  const unsigned int d4_c ,
		  const unsigned int d5_c ,
		  const unsigned int d6_c ,
		  const unsigned int d7_c ,
		  const unsigned int d8_c) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate (d0_c , d1_c , d2_c , d3_c , d4_c , d5_c , d6_c , d7_c , d8_c);
  }

  array (const class array<C> &X) : d0 (0) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (0) , table (NULL)
  {
    allocate_fill (X);
  }

  explicit array (const unsigned int N , const C T[]) : d0 (N) , d1 (0) , d2 (0) , d3 (0) , d4 (0) , d5 (0) , d6 (0) , d7 (0) , d8 (0) , dimension_table (N) , table (NULL)
  {
    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] = T[i];
  }

  ~array ()
  {
    delete [] table;
  }

  void allocate (const unsigned int d0_c)
  {
    // if (dimension_table > 0) {
    //   //std::cerr<<"array cannot be allocated twice in array::allocate (1)"<<std::endl;
    //   exit(0);
    // }

    d0 = d0_c , dimension_table = d0;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (2)");

    d0 = d0_c , d1 = d1_c , dimension_table = d0 * d1;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (3)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , dimension_table = d0 * d1 * d2;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (4)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , dimension_table = d0 * d1 * d2 * d3;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c ,
		 const unsigned int d4_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (5)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , d4 = d4_c , dimension_table = d0 * d1 * d2 * d3 * d4;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c ,
		 const unsigned int d4_c ,
		 const unsigned int d5_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (6)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , d4 = d4_c , d5 = d5_c , dimension_table = d0 * d1 * d2 * d3 * d4 * d5;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c ,
		 const unsigned int d4_c ,
		 const unsigned int d5_c ,
		 const unsigned int d6_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (7)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , d4 = d4_c , d5 = d5_c , d6 = d6_c , dimension_table = d0 * d1 * d2 * d3 * d4 * d5 * d6;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c ,
		 const unsigned int d4_c ,
		 const unsigned int d5_c ,
		 const unsigned int d6_c ,
		 const unsigned int d7_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (8)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , d4 = d4_c , d5 = d5_c , d6 = d6_c , d7 = d7_c , dimension_table = d0 * d1 * d2 * d3 * d4 * d5 * d6 * d7;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate (
		 const unsigned int d0_c ,
		 const unsigned int d1_c ,
		 const unsigned int d2_c ,
		 const unsigned int d3_c ,
		 const unsigned int d4_c ,
		 const unsigned int d5_c ,
		 const unsigned int d6_c ,
		 const unsigned int d7_c ,
		 const unsigned int d8_c)
  {
    if (dimension_table > 0) error_message_print_abort ("array cannot be allocated twice in array::allocate (9)");

    d0 = d0_c , d1 = d1_c , d2 = d2_c , d3 = d3_c , d4 = d4_c , d5 = d5_c , d6 = d6_c , d7 = d7_c , d8 = d8_c , dimension_table = d0 * d1 * d2 * d3 * d4 * d5 * d6 * d7 * d8;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);
  }

  void allocate_fill (const class array<C> &X)
  {
    d0 = X.d0 , d1 = X.d1 , d2 = X.d2 , d3 = X.d3 , d4 = X.d4 , d5 = X.d5 , d6 = X.d6 , d7 = X.d7 , d8 = X.d8 , dimension_table = X.dimension_table;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] = X.table[i];
  }

  // copy allocator for objects necessitating an explicit allocation

  void allocate_fill_object_elements (const class array<C> &X)
  {
    d0 = X.d0 , d1 = X.d1 , d2 = X.d2 , d3 = X.d3 , d4 = X.d4 , d5 = X.d5 , d6 = X.d6 , d7 = X.d7 , d8 = X.d8 , dimension_table = X.dimension_table;

    table = (dimension_table > 0) ? (new C [dimension_table]) : (NULL);

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i].allocate_fill (X.table[i]);
  }

  void deallocate ()
  {
    delete [] table;

    d0 = d1 = d2 = d3 = d4 = d5 = d6 = d7 = d8 = dimension_table = 0;

    table = NULL;
  }

  // deallocator deallocating only objects elements , not the array itself

  void deallocate_object_elements ()
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i].deallocate ();
  }

  bool is_it_filled () const
  {
    return (dimension_table > 0);
  }

  void operator = (const class array<C> &X)
  {
    const C *const X_table = X.table;

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] = X_table[i];
  }

  void operator = (const C &x)
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] = x;
  }

  void operator += (const C &x)
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] += x;
  }

  void operator -= (const C &x)
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] -= x;
  }

  void operator *= (const C &x)
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] *= x;
  }

  void operator /= (const C &x)
  {
    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] /= x;
  }

  void operator += (const class array<C> &A)
  {
    const C *const A_table = A.table;

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] += A_table[i];
  }

  void operator -= (const class array<C> &A)
  {
    const C *const A_table = A.table;

    for (unsigned int i = 0 ; i < dimension_table ; i++) table[i] -= A_table[i];
  }

  inline C & operator [] (const unsigned int i) const
  {
    return table[i];
  }

  C & operator () (const unsigned int i0) const
  {
    return table[i0];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1) const
  {
    return table[i1 + d1 * i0];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2) const
  {
    return table[i2 + d2 * (i1 + d1 * i0)];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3) const
  {
    return table[i3 + d3 * (i2 + d2 * (i1 + d1 * i0))];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3 ,
		   const unsigned int i4) const
  {
    return table[i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)))];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3 ,
		   const unsigned int i4 ,
		   const unsigned int i5) const
  {
    return table[i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0))))];
  }

  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3 ,
		   const unsigned int i4 ,
		   const unsigned int i5 ,
		   const unsigned int i6) const
  {
    return table[i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)))))];
  }


  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3 ,
		   const unsigned int i4 ,
		   const unsigned int i5 ,
		   const unsigned int i6 ,
		   const unsigned int i7) const
  {
    return table[i7 + d7 * (i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0))))))];
  }


  C & operator () (
		   const unsigned int i0 ,
		   const unsigned int i1 ,
		   const unsigned int i2 ,
		   const unsigned int i3 ,
		   const unsigned int i4 ,
		   const unsigned int i5 ,
		   const unsigned int i6 ,
		   const unsigned int i7 ,
		   const unsigned int i8) const
  {
    return table[i8 + d8*(i7 + d7 * (i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)))))))];
  }

  // One-dimensional index calculated from its N-dimensional indices
  // To be used with operator [] for example.

  unsigned int index_determine (const unsigned int i0) const
  {
    return i0;
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1) const
  {
    return (i1 + d1 * i0);
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2) const
  {
    return (i2 + d2 * (i1 + d1 * i0));
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3) const
  {
    return (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)));
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3 ,
				const unsigned int i4) const
  {
    return (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0))));
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3 ,
				const unsigned int i4 ,
				const unsigned int i5) const
  {
    return (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)))));
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3 ,
				const unsigned int i4 ,
				const unsigned int i5 ,
				const unsigned int i6) const
  {
    return (i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0))))));
  }


  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3 ,
				const unsigned int i4 ,
				const unsigned int i5 ,
				const unsigned int i6 ,
				const unsigned int i7) const
  {
    return (i7 + d7 * (i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0)))))));
  }

  unsigned int index_determine (
				const unsigned int i0 ,
				const unsigned int i1 ,
				const unsigned int i2 ,
				const unsigned int i3 ,
				const unsigned int i4 ,
				const unsigned int i5 ,
				const unsigned int i6 ,
				const unsigned int i7 ,
				const unsigned int i8) const
  {
    return (i8 + d8*(i7 + d7 * (i6 + d6 * (i5 + d5 * (i4 + d4 * (i3 + d3 * (i2 + d2 * (i1 + d1 * i0))))))));
  }

  void read_disk (const string &file_name)
  {
    basic_read_disk<C> (file_name , dimension_table , table);
  }

  void copy_disk (const string &file_name) const
  {
    basic_copy_disk<C> (file_name , dimension_table , table);
  }

  C sum () const
  {
    C sum_table = 0.0;

    for (unsigned int i = 0 ; i < dimension_table ; i++) sum_table += table[i];

    return sum_table;
  }

  C product () const
  {
    C product_table = 1.0;

    for (unsigned int i = 0 ; i < dimension_table ; i++) product_table *= table[i];

    return product_table;
  }

  unsigned int dimension (const unsigned int i) const
  {
    if (i == 0) return d0;
    if (i == 1) return d1;
    if (i == 2) return d2;
    if (i == 3) return d3;
    if (i == 4) return d4;
    if (i == 5) return d5;
    if (i == 6) return d6;
    if (i == 7) return d7;
    if (i == 8) return d8;

    return 0;
  }

  unsigned int dimension_total () const
  {
    return dimension_table;
  }

  C max () const;
  C min () const;

  void quick_sort (int low , int high);
  void quick_sort_by_real_part (int low , int high);
  void quick_sort_by_abs (int low , int high);

  void reverse_order ();

  void assign (const class array<C> &T);

  unsigned int closest_value_index_determine (const C &x) const;
  bool is_it_equally_spaced () const;

  unsigned int first_index_determine_for_MPI (const unsigned int group_processes_number , const unsigned int process) const;
  unsigned int last_index_determine_for_MPI  (const unsigned int group_processes_number , const unsigned int process) const;

  unsigned int active_process_determine_for_MPI (const unsigned int group_processes_number , const unsigned int index) const;

#ifdef UseMPI
  void MPI_Send  (const unsigned int Recv_process , const int tag , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Send  (const unsigned int Recv_process , const int tag , const MPI_Comm MPI_C);

  void MPI_Recv  (const unsigned int Send_process , const int tag , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Recv  (const unsigned int Send_process , const int tag , const MPI_Comm MPI_C);

  void MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Comm MPI_C);

  void MPI_Bcast (const unsigned int Send_process , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Bcast (const unsigned int Send_process , const MPI_Comm MPI_C);

  void MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Comm MPI_C);

  void MPI_Reduce    (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Reduce    (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Comm MPI_C);

  void MPI_Allreduce (MPI_Op op , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C);
  void MPI_Allreduce (MPI_Op op , const MPI_Comm MPI_C);
#endif

private:

  unsigned int d0 , d1 , d2 , d3 , d4 , d5 , d6 , d7 , d8 , dimension_table; // dimensions of the array.

  C * table; // components of the array.
};

// Memory stored in the array in Mb
// --------------------------------

template <class C>
double used_memory_calc (const class array<C> &A)
{
  double used_memory = sizeof (class array<C>)/1000000.0;

  for (unsigned int i = 0 ; i < A.dimension_total () ; i++) used_memory += used_memory_calc (A[i]);

  return used_memory;
}


// Maximal value of array
// ----------------------
template <class C>
C array<C>::max () const
{
  if (dimension_table == 0) error_message_print_abort ("max function not defined for arrays without elements");

  C table_max = table[0];

  for (unsigned int i = 1 ; i < dimension_table ; i++)
    {
      const C table_i = table[i];

      if (table_max < table_i) table_max = table_i;
    }

  return table_max;
}

// Minimal value of array
// ----------------------
template <class C>
C array<C>::min () const
{
  if (dimension_table == 0) error_message_print_abort ("min function not defined for arrays without elements");

  C table_min = table[0];

  for (unsigned int i = 1 ; i < dimension_table ; i++)
    {
      const C table_i = table[i];

      if (table_min > table_i) table_min = table_i;
    }

  return table_min;
}


// Quick sort of array (general)
// -----------------------------
template<class C>
void array<C>::quick_sort (int low , int high)
{
  if (dimension_table == 0) return;

  if (d1 > 0) error_message_print_abort ("Only one-dimensional arrays can be sorted (general).");

  ::quick_sort<C> (table , low , high);
}



// Quick sort of array (by real parts)
// -----------------------------------
template<class C>
void array<C>::quick_sort_by_real_part (int low , int high)
{
  if (dimension_table == 0) return;

  if (d1 > 0) error_message_print_abort ("Only one-dimensional arrays can be sorted (real part).");

  ::quick_sort_by_real_part<C> (table , low , high);
}



// Quick sort of array (by moduli)
// -------------------------------
template<class C>
void array<C>::quick_sort_by_abs (int low , int high)
{
  if (dimension_table == 0) return;

  if (d1 > 0) error_message_print_abort ("Only one-dimensional arrays can be sorted (abs part).");

  ::quick_sort_by_abs<C> (table , low , high);
}





// Reverses the order of the array
// -------------------------------
// {a , b , c , d} becomes {d , c , b , a} for example
// Valid for 1D table only.

template<class C>
void array<C>::reverse_order ()
{
  if ((dimension_table != d0) || (d1 != 0) || (d2 != 0) || (d3 != 0) || (d4 != 0) || (d5 != 0) || (d6 != 0) || (d7 != 0) || (d8 != 0))
    error_message_print_abort ("reverse_order in array class only for 1D arrays.");

  if (dimension_table == 0) return;

  const unsigned int half_dimension_table = dimension_table/2;

  const unsigned int dimension_table_minus_one = dimension_table - 1;

  for (unsigned int i = 0 ; i < half_dimension_table ; i++) swap<C> (table[i] , table[dimension_table_minus_one - i]);
}


// Assignment
// ----------
// T is assigned to A if T.dimension > A.dimension , that is A(i) = T(i) when available while remaining components of T are ignored.
// T is copied to A if T.dimension < A.dimension , remaining components of A being ignored.
// 1D tables only.

template <typename C>
void array<C>::assign (const class array<C> &T)
{
  if ((dimension_table != d0) || (d1 != 0) || (d2 != 0) || (d3 != 0) || (d4 != 0) || (d5 != 0) || (d6 != 0) || (d7 != 0) || (d8 != 0))
    error_message_print_abort ("assign in array class only for 1D arrays.");

  class array<C> &A = *this;

  const unsigned int smallest_dimension_table = ::min (dimension_table , T.dimension_table);

  for (unsigned int i = 0 ; i < smallest_dimension_table ; i++) A(i) = T(i);
}


// Index of the component closest to a given number
// ------------------------------------------------
template<typename C>
unsigned int array<C>::closest_value_index_determine (const C &x) const
{
  if ((dimension_table == 0) || (d1 > 0)) error_message_print_abort ("closest_value_index_determine for non empty one-dimensional arrays only.");

  double difference = abs (table[0] - x);

  unsigned int closest_value_index = 0;

  for (unsigned int i = 1 ; i < dimension_table ; i++)
    {
      const double difference_i = abs (table[i] - x);

      if (difference_i < difference)
	{
	  closest_value_index = i;

	  difference = difference_i;
	}
    }

  return closest_value_index;
}

// determine if an array is equally spaced
// ---------------------------------------
// true if T[i] = T[0] + i.step , false if not

template <typename SCALAR_TYPE>
bool array<SCALAR_TYPE>::is_it_equally_spaced () const
{
  if (dimension_table <= 1) return true;

  if (d1 > 0) error_message_print_abort ("closest_value_index_determine for one-dimensional arrays only.");

  const SCALAR_TYPE step = table[1] - table[0];

  const unsigned int dimension_table_minus_one = dimension_table - 1;

  for (unsigned int i = 1 ; i < dimension_table_minus_one ; i++)
    {
      if (abs (table[i + 1] - table[i] - step) > precision) return false;
    }

  return true;
}






// Unary "+" overloading.
// ----------------------
template <typename C>
class array<C> operator + (const class array<C> &A)
{
  return A;
}



// Unary "-" overloading.
// ----------------------
template <typename C>
class array<C> operator - (const class array<C> &A)
{
  const unsigned int dimension = A.dimension_total ();

  class array<C> M(dimension);

  for (unsigned int i = 0 ; i < dimension ; i++) M[i] = -A[i];

  return M;
}







// "+" overloading.
// ----------------
template <typename C>
class array<C> operator + (const class array<C> &A , const class array<C> &B)
{
  class array<C> M = A;

  M += B;

  return M;
}

template <typename C>
class array<C> operator + (const class array<C> &A , const C &x)
{
  class array<C> M = A;

  M += x;

  return M;
}


template <typename C>
class array<C> operator + (const C &x , const class array<C> &A)
{
  class array<C> M = A;

  M += x;

  return M;
}




// "-" overloading.
// ----------------
template <typename C>
class array<C> operator - (const class array<C> &A , const class array<C> &B)
{
  class array<C> M = A;

  M -= B;

  return M;
}

template <typename C>
class array<C> operator - (const class array<C> &A , const C &x)
{
  class array<C> M = A;

  M -= x;

  return M;
}


template <typename C>
class array<C> operator - (const C &x , const class array<C> &A)
{
  class array<C> M = -A;

  M += x;

  return M;
}



// "*" overloading : A times x
//----------------------------
template <typename C>
class array<C> operator * (const class array<C> &A , const C &x)
{
  class array<C> M = A;

  M *= x;

  return M;
}



// "*" overloading : x times A
// ---------------------------
template <typename C>
class array<C> operator * (const C &x , const class array<C> &A)
{
  class array<C> M = A;

  M *= x;

  return M;
}





// "/" overloading : A over x
// --------------------------
template <typename C>
class array<C> operator / (const class array<C> &A , const C &x)
{
  class array<C> M = A;

  M /= x;

  return M;
}



// "==" overloading
// ----------------
template <typename C>
bool operator == (const class array<C> &A , const class array<C> &B)
{
  const unsigned int dimension_table = A.dimension_total ();

  for (unsigned int i = 0 ; i < dimension_table ; i++)
    {
      if (A[i] != B[i]) return false;
    }

  return true;
}

// "!=" overloading
// ----------------
template <typename C>
bool operator != (const class array<C> &A , const class array<C> &B)
{
  return (!(A == B));
}


// Standard MPI routines (see MPI_helper)
// --------------------------------------


template <class C>
unsigned int array<C>::first_index_determine_for_MPI (const unsigned int group_processes_number , const unsigned int process) const
{
  return basic_first_index_determine_for_MPI (dimension_table , group_processes_number , process);
}

template <class C>
unsigned int array<C>::last_index_determine_for_MPI (const unsigned int group_processes_number , const unsigned int process) const
{
  return basic_last_index_determine_for_MPI (dimension_table , group_processes_number , process);
}

template <class C>
unsigned int array<C>::active_process_determine_for_MPI (const unsigned int group_processes_number , const unsigned int index) const
{
  return basic_active_process_determine_for_MPI (dimension_table , group_processes_number , index);
}



#ifdef UseMPI

template <class C>
void array<C>::MPI_Send (const unsigned int Recv_process , const int tag , 	const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Send<C> (dimension_table , table , Recv_process , tag , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Send (const unsigned int Recv_process , const int tag , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Send<C> (dimension_table , table , Recv_process , tag , MPI_C);
}

template <class C>
void array<C>::MPI_Recv (const unsigned int Send_process , const int tag , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Recv<C> (dimension_table , table , Send_process , tag , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Recv (const unsigned int Send_process , const int tag , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Recv<C> (dimension_table , table , Send_process , tag , MPI_C);
}

template <class C>
void array<C>::MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Sendrecv_replace<C> (dimension_table , table , Send_process , Recv_process , Send_tag , Recv_tag ,  MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Sendrecv_replace<C> (dimension_table , table , Send_process , Recv_process , Send_tag , Recv_tag ,  MPI_C);
}

template <class C>
void array<C>::MPI_Bcast (const unsigned int Send_process , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Bcast<C> (dimension_table , table ,  Send_process , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Bcast (const unsigned int Send_process , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Bcast<C> (dimension_table , table , Send_process , MPI_C);
}

template <class C>
void array<C>::MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Allgatherv<C> (dimension_table , table , group_processes_number , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Allgatherv<C> (dimension_table , table , group_processes_number , MPI_C);
}

template <class C>
void array<C>::MPI_Reduce (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Reduce<C> (dimension_table , table , op , Recv_process , process , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Reduce (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Reduce<C> (dimension_table , table , op , Recv_process , process , MPI_C);
}

template <class C>
void array<C>::MPI_Allreduce (MPI_Op op , const MPI_Datatype MPI_T_type , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Allreduce<C> (dimension_table , table , op , MPI_T_type , MPI_C);
}

template <class C>
void array<C>::MPI_Allreduce (MPI_Op op , const MPI_Comm MPI_C)
{
  if (dimension_table == 0) return;

  MPI_helper::Allreduce<C> (dimension_table , table , op , MPI_C);
}

#endif





//  array print overload (1D)
// --------------------------
template <typename C>
ostream & operator << (ostream &os , const class array<C> &T)
{
  const unsigned int N = T.dimension_total ();

  for (unsigned int i = 0 ; i < N ; i++) os << T[i] << " ";

  return os;
}



// Class describing a vector type.
// -------------------------------
// SCALAR_TYPE must be double or complex.
// One cannot mix different scalar types in routines.
// Use convert or complex_vector to go from one scalar type to an another.
// List is a standard allocated pointer.
//
// One can test if a vector is allocate with is_it_filled
//
// Standard operator overloadings of the form = , + , - , * , / and -= , += , *= , /= can be used.
// Comparison operators of form == , != can be used.
//
// One can calculate the ||oo norm and the 2-norm with infinite_norm and dot_product_norm.
//
// One can normalize the vector with normalization
//
// One can assign the content of a vector in a vector of different dimension with assign
//
// One can generate random vectors with random_vector (arbitrary seed) and real_random_vector (arbitrary seed , imaginary parts are zero) ,
// pseudo_random_vector (seed is the dimension of the vector) , pseudo_real_random_vector (seed is the dimension of the vector , imaginary parts are zero) ,
//
// One can calculate the complex conjugate of a vector with conjugate
//
// One can fix its phase with good_phase (the largest component real or imaginary part , whatever is larger , gets positive sign)
//
// One can test if a vector is finite with isfinite
//
// One can test if a vector is real with is_it_real
//
// One can reverse the order of elements with reverse_order
//
// One can find the index of a value closest to a given value with closest_value_index_determine
//
// Standard MPI routines are implemented (see MPI_helper)
//
// One can copy and read a vector from disk with copy_disk (output.dat) and read_disk (input.dat) for example. Files are not formatted.

template <typename SCALAR_TYPE>
class vector_class
{
public :

  explicit vector_class () : dimension (0) , List (NULL) {}

  explicit vector_class (const unsigned int dimension_c) : dimension (dimension_c) , List (NULL)
  {
    List = (dimension > 0) ? (new SCALAR_TYPE [dimension]) : (NULL);
  }

  vector_class (const class vector_class &X) : dimension (X.dimension) , List (NULL)
  {
    List = (dimension > 0) ? (new SCALAR_TYPE [dimension]) : (NULL);

    for (unsigned int i = 0 ; i < dimension ; i++) List[i] = X.List[i];
  }

  ~vector_class ()
  {
    delete [] List;
  }

  unsigned int get_dimension () const
  {
    return dimension;
  }

  void allocate (const unsigned int dimension_c)
  {
    if (dimension > 0) error_message_print_abort ("vector_class cannot be allocated twice in vector_class::allocate");

    dimension = dimension_c;

    delete [] List;

    List = (dimension > 0) ? (new SCALAR_TYPE [dimension]) : (NULL);
  }

  void allocate_fill (const class vector_class &X)
  {
    dimension = X.dimension;

    delete [] List;

    List = (dimension > 0) ? (new SCALAR_TYPE [dimension]) : (NULL);

    for (unsigned int i = 0 ; i < dimension ; i++) List[i] = X.List[i];
  }

  void deallocate ()
  {
    dimension = 0;

    delete [] List;

    List = NULL;
  }

  bool is_it_filled () const
  {
    return (dimension > 0);
  }

  void read_disk (const string &file_name)
  {
    basic_read_disk<SCALAR_TYPE> (file_name , dimension , List);
  }

  void copy_disk (const string &file_name) const
  {
    basic_copy_disk<SCALAR_TYPE> (file_name , dimension , List);
  }

  SCALAR_TYPE sum () const
  {
    SCALAR_TYPE sum_table = 0.0;

    for (unsigned int i = 0 ; i < dimension ; i++) sum_table += List[i];

    return sum_table;
  }

  SCALAR_TYPE product () const
  {
    SCALAR_TYPE product_table = 1.0;

    for (unsigned int i = 0 ; i < dimension ; i++) product_table *= List[i];

    return product_table;
  }

  SCALAR_TYPE max () const;
  SCALAR_TYPE min () const;

  class vector_class& operator += (const class vector_class &);
  class vector_class& operator -= (const class vector_class &);

  class vector_class& operator += (const SCALAR_TYPE &);
  class vector_class& operator -= (const SCALAR_TYPE &);
  class vector_class& operator *= (const SCALAR_TYPE &);
  class vector_class& operator /= (const SCALAR_TYPE &);

  const class vector_class& operator = (const SCALAR_TYPE &);
  const class vector_class& operator = (const vector_class &);

  SCALAR_TYPE & operator() (const unsigned int) const;
  SCALAR_TYPE & operator[] (const unsigned int) const;

  double infinite_norm () const;
  SCALAR_TYPE dot_product_norm () const;
  void normalization ();
  void assign (const class vector_class<SCALAR_TYPE> &V);
  void random_vector ();
  void real_random_vector ();
  void pseudo_random_vector ();
  void pseudo_real_random_vector ();

  void conjugate ();
  void good_phase ();
  bool isfinite () const;
  bool is_it_real () const;

  void quick_sort (int low , int high);
  void quick_sort_by_real_part (int low , int high);
  void quick_sort_by_abs (int low , int high);

  void reverse_order ();

  unsigned int closest_value_index_determine (const SCALAR_TYPE &x) const;

#ifdef UseMPI
  void MPI_Send  (const unsigned int Recv_process , const int tag , const MPI_Comm MPI_C);
  void MPI_Recv  (const unsigned int Send_process , const int tag , const MPI_Comm MPI_C);
  void MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Comm MPI_C);
  void MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Comm MPI_C);
  void MPI_Bcast (const unsigned int Send_process , const MPI_Comm MPI_C);
  void MPI_Reduce (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Comm MPI_C);
  void MPI_Allreduce (MPI_Op op , const MPI_Comm MPI_C);
#endif

private:

  unsigned int dimension;  // dimension of the vector.

  SCALAR_TYPE *List;       // components of the vector.
};






// "()" overloading.
// -----------------
// V(i) = V.List[i]

template <typename SCALAR_TYPE>
SCALAR_TYPE & vector_class<SCALAR_TYPE>::operator() (const unsigned int i) const
{
  return List[i];
}

template <typename SCALAR_TYPE>
SCALAR_TYPE & vector_class<SCALAR_TYPE>::operator[] (const unsigned int i) const
{
  return List[i];
}




// "=" overloading.
// ----------------
template <typename SCALAR_TYPE>
const class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator = (const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) = V(i);

  return W;
}


// "=" overloading.
// ----------------
// All components are equal to the same value

template <typename SCALAR_TYPE>
const class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator = (const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) = x;

  return W;
}


// "+=" overloading.
// ----------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator += (const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) += V(i);

  return W;
}


// "-=" overloading.
// ----------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator -= (const class vector_class<SCALAR_TYPE> & V)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) -= V(i);

  return W;
}




// "+=" overloading : vector plus SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator += (const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) += x;

  return W;
}


// "-=" overloading : vector minus SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator -= (const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) -= x;

  return W;
}

// "*=" overloading : vector times SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator *= (const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) *= x;

  return W;
}


// "/=" overloading : vector over SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> & vector_class<SCALAR_TYPE>::operator /= (const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) /= x;

  return W;
}







// "+" overloading.
// ----------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator + (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const class vector_class<SCALAR_TYPE> &W)
{
  class vector_class<SCALAR_TYPE> S = V;

  S += W;

  return S;
}

template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator + (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> S = V;

  S += x;

  return S;
}

template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator + (
					    const SCALAR_TYPE &x ,
					    const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> S = V;

  S += x;

  return S;
}



// "-" overloading.
// ----------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator - (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const class vector_class<SCALAR_TYPE> &W)
{
  class vector_class<SCALAR_TYPE> S = V;

  S -= W;

  return S;
}

template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator - (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> S = V;

  S -= x;

  return S;
}


template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator - (
					    const SCALAR_TYPE &x ,
					    const class vector_class<SCALAR_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  class vector_class<SCALAR_TYPE> S(N);

  S = x;

  S -= V;

  return S;
}





// "*=" overloading : vector times SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator * (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> W = V;

  W *= x;

  return W;
}




// "*" overloading : vector times SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator * (
					    const SCALAR_TYPE &x ,
					    const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> W = V;

  W *= x;

  return W;
}





// "/" overloading : vector over SCALAR_TYPE
// ------------------------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator / (
					    const class vector_class<SCALAR_TYPE> &V ,
					    const SCALAR_TYPE &x)
{
  class vector_class<SCALAR_TYPE> W = V;

  W /= x;

  return W;
}






// Infinite norm.
// --------------
// It is the largest component |V(i)|oo , with i in [0:V.dimension-1].

template <typename SCALAR_TYPE>
double vector_class<SCALAR_TYPE>::infinite_norm () const
{
  const class vector_class<SCALAR_TYPE> &W = *this;

  double inf_norm_value = 0.0;

  for (unsigned int i = 0 ; i < dimension ; i++)
    {
      const double inf_norm_Wi = inf_norm (W(i));

      if (!finite (inf_norm_Wi)) return inf_norm_Wi;

      if (inf_norm_value < inf_norm_Wi) inf_norm_value = inf_norm_Wi;
    }

  return inf_norm_value;
}





//  norm 2 : not the hermitian one for complex , but the Berggren one.
// ------------------------------------------------------------------
// It is sqrt[\sum_{i=0}^{V.dimension-1} V(i)^2].

template <typename SCALAR_TYPE>
SCALAR_TYPE vector_class<SCALAR_TYPE>::dot_product_norm () const
{
  const class vector_class<SCALAR_TYPE> &W = *this;

  return sqrt (W*W);
}




// Unary +
// -------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator + (const class vector_class<SCALAR_TYPE> &V)
{
  return V;
}


// Unary -
// -------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> operator - (const class vector_class<SCALAR_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  class vector_class<SCALAR_TYPE> W(N);

  for (unsigned int i = 0 ; i < N ; i++) W(i) = -V(i);

  return W;
}




// Dot product : not the hermitian one for complex , but the Berggren one.
// ----------------------------------------------------------------------
// It is \sum_{i=0}^{V.dimension-1} V(i)*W(i).


template <typename SCALAR_TYPE>
SCALAR_TYPE operator * (
			const class vector_class<SCALAR_TYPE> &V ,
			const class vector_class<SCALAR_TYPE> &W)
{
  const unsigned int N = V.get_dimension ();

  SCALAR_TYPE overlap = 0.0;

  for (unsigned int i = 0 ; i < N ; i++) overlap += V(i)*W(i);

  return overlap;
}





//  normalization with the Berggren norm.
// --------------------------------------
template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::normalization ()
{
  class vector_class<SCALAR_TYPE> &W = *this;

  W /= dot_product_norm ();
}






// Assignment
// ----------
// V is assigned to W if V.dimension > W.dimension , that is W(i) = V(i) when available while remaining components of V are ignored.
// V is copied to W if V.dimension < W.dimension , remaining components of W being gnored.

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::assign (const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> &W = *this;

  const unsigned int smallest_dimension = ::min (dimension , V.get_dimension ());

  for (unsigned int i = 0 ; i < smallest_dimension ; i++) W(i) = V(i);
}




//  V randomly initialized.
// ------------------------

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::real_random_vector ()
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) = random_number<double> () - 0.5;
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::random_vector ()
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) = random_number<SCALAR_TYPE> () - 0.5;
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::pseudo_real_random_vector ()
{
  seed (0);

  real_random_vector ();

  seed ();
}


template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::pseudo_random_vector ()
{
  seed (0);

  random_vector ();

  seed ();
}





#ifdef UseMPI

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Send (const unsigned int Recv_process , const int tag , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Send<SCALAR_TYPE> (dimension , List , Recv_process , tag , MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Recv (const unsigned int Send_process , const int tag , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Recv<SCALAR_TYPE> (dimension , List , Send_process , tag , MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Sendrecv_replace (const unsigned int Send_process , const unsigned int Recv_process , const int Send_tag , const int Recv_tag , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Sendrecv_replace<SCALAR_TYPE> (dimension , List , Send_process , Recv_process , Send_tag , Recv_tag ,  MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Allgatherv (const unsigned int group_processes_number , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Allgatherv<SCALAR_TYPE> (dimension , List , group_processes_number , MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Bcast (const unsigned int Send_process , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Bcast<SCALAR_TYPE> (dimension , List , Send_process , MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Reduce (MPI_Op op , const unsigned int Recv_process , const unsigned int process , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Reduce<SCALAR_TYPE> (dimension , List , op , Recv_process , process , MPI_C);
}

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::MPI_Allreduce (MPI_Op op , const MPI_Comm MPI_C)
{
  if (dimension == 0) return;

  MPI_helper::Allreduce<SCALAR_TYPE> (dimension , List , op , MPI_C);
}

#endif









//  V << print overload.
// ---------------------
template <typename SCALAR_TYPE>
ostream & operator << (ostream &os , const class vector_class<SCALAR_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  for (unsigned int i = 0 ; i < N ; i++) os << V(i) << " ";

  return os;
}




//  V >> get overload.
// -------------------
template <typename SCALAR_TYPE>
istream & operator >> (istream &is , class vector_class<SCALAR_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  for (unsigned int i = 0 ; i < N ; i++) is >> V(i);

  return is;
}

// Vector conjugate
// ----------------
template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::conjugate ()
{
  class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++) W(i) = conj_dc (W(i));
}


// Phase of vector fixed
// ---------------------
template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::good_phase ()
{
  class vector_class<SCALAR_TYPE> &W = *this;

  double inf_norm_value = 0.0;

  unsigned int i_max = 0;

  for (unsigned int i = 0 ; i < dimension ; i++)
    {
      const double try_inf_norm_value = inf_norm (W(i));

      if (inf_norm_value < try_inf_norm_value)
	{
	  inf_norm_value = try_inf_norm_value;

	  i_max = i;
	}
    }

  const double re_max = real_dc (W(i_max));
  const double im_max = imag_dc (W(i_max));

  const int phase = (abs (re_max) > abs (im_max)) ? (SIGN (re_max)) : (SIGN (im_max));

  if (phase == -1)
    {
      for (unsigned int i = 0 ; i < dimension ; i++) W(i) = -W(i);
    }
}

// Test if vector is finite
// ------------------------
template <typename SCALAR_TYPE>
bool vector_class<SCALAR_TYPE>::isfinite () const
{
  const class vector_class<SCALAR_TYPE> &W = *this;

  for (unsigned int i = 0 ; i < dimension ; i++)
    {
      if (!finite (W(i))) return false;
    }

  return true;
}



// Check if V is a real vector
// ---------------------------
template <typename SCALAR_TYPE>
bool vector_class<SCALAR_TYPE>::is_it_real () const
{
  const class vector_class<SCALAR_TYPE> &V = *this;

  for (unsigned int i = 0 ; i < dimension ; i++)
    {
      if (imag_dc (V(i)) != 0) return false;
    }

  return true;
}



// Reverses the order of the vector
// --------------------------------
// {a , b , c , d} becomes {d , c , b , a} for example

template <typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::reverse_order ()
{
  if (dimension == 0) return;

  const unsigned int half_dimension = dimension/2;

  const unsigned int dimension_minus_one = dimension - 1;

  class vector_class<SCALAR_TYPE> &V = *this;

  for (unsigned int i = 0 ; i < half_dimension ; i++) swap<SCALAR_TYPE> (V(i) , V(dimension_minus_one - i));
}



// Quick sort of array (general)
// -----------------------------
template<typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::quick_sort (int low , int high)
{
  if (dimension == 0) return;

  ::quick_sort<SCALAR_TYPE> (List , low , high);
}



// Quick sort of array (by real parts)
// -----------------------------------
template<typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::quick_sort_by_real_part (int low , int high)
{
  if (dimension == 0) return;

  ::quick_sort_by_real_part<SCALAR_TYPE> (List , low , high);
}



// Quick sort of array (by moduli)
// -------------------------------
template<typename SCALAR_TYPE>
void vector_class<SCALAR_TYPE>::quick_sort_by_abs (int low , int high)
{
  if (dimension == 0) return;

  ::quick_sort_by_abs<SCALAR_TYPE> (List , low , high);
}








// Index of value closest to an initial given value
// ------------------------------------------------
template <typename SCALAR_TYPE>
unsigned int vector_class<SCALAR_TYPE>::closest_value_index_determine (const SCALAR_TYPE &x) const
{
  if (dimension == 0) error_message_print_abort ("closest_value_index_determine for non empty vectors only.");

  const class vector_class<SCALAR_TYPE> &V = *this;

  double difference = abs (V(0) - x);

  unsigned int closest_value_index = 0;

  for (unsigned int i = 1 ; i < dimension ; i++)
    {
      const double difference_i = abs (V(i) - x);

      if (difference_i < difference)
	{
	  closest_value_index = i;

	  difference = difference_i;
	}
    }

  return closest_value_index;
}



// Vector real part returned
// -------------------------
template <typename SCALAR_TYPE , typename COMPLEX_TYPE>
class vector_class<SCALAR_TYPE> real (const class vector_class<COMPLEX_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  class vector_class<SCALAR_TYPE> W(N);

  for (unsigned int i = 0 ; i < N ; i++) W(i) = real_dc (V(i));

  return W;
}

// Vector imaginary part returned
// ------------------------------
template <typename SCALAR_TYPE , typename COMPLEX_TYPE>
class vector_class<SCALAR_TYPE> imag (const class vector_class<COMPLEX_TYPE> &V)
{
  const unsigned int N = V.get_dimension ();

  class vector_class<SCALAR_TYPE> W(N);

  for (unsigned int i = 0 ; i < N ; i++) W(i) = imag_dc (V(i));

  return W;
}

// Vector conjugate returned
// -------------------------
template <typename SCALAR_TYPE>
class vector_class<SCALAR_TYPE> conj (const class vector_class<SCALAR_TYPE> &V)
{
  class vector_class<SCALAR_TYPE> W = V;

  W.conjugate ();

  return W;
}


// Conversion of one scalar type to another
// ----------------------------------------

template <typename IN_SCALAR_TYPE , typename OUT_SCALAR_TYPE>
class vector_class<OUT_SCALAR_TYPE> convert (const class vector_class<IN_SCALAR_TYPE> &V_in)
{
  const unsigned int N = V_in.get_dimension ();

  class vector_class<OUT_SCALAR_TYPE> V_out(N);

  for (unsigned int i = 0 ; i < N ; i++) V_out(i) = V_in(i);

  return V_out;
}




// Complex vector returned from real and/or imaginary parts
// --------------------------------------------------------

template <typename SCALAR_TYPE , typename COMPLEX_TYPE>
class vector_class<COMPLEX_TYPE> complex_vector (const class vector_class<SCALAR_TYPE> &Re_V)
{
  const unsigned int N = Re_V.get_dimension ();

  class vector_class<COMPLEX_TYPE> V(N);

  for (unsigned int i = 0 ; i < N ; i++) V(i) = Re_V(i);

  return V;
}

template <typename SCALAR_TYPE , typename COMPLEX_TYPE>
class vector_class<COMPLEX_TYPE> complex_vector (
						 const class vector_class<SCALAR_TYPE> &Re_V ,
						 const class vector_class<SCALAR_TYPE> &Im_V)
{
  const unsigned int N = Re_V.get_dimension ();

  class vector_class<COMPLEX_TYPE> V(N);

  for (unsigned int i = 0 ; i < N ; i++) V(i) = COMPLEX_TYPE (Re_V(i) , Im_V(i));

  return V;
}





template <typename SCALAR_TYPE>
bool operator == (const class vector_class<SCALAR_TYPE> &V , const class vector_class<SCALAR_TYPE> &W)
{
  const unsigned int dimension = V.get_dimension ();

  for (unsigned int i = 0 ; i < dimension ; i++)
    {
      if (V(i) != W(i)) return false;
    }

  return true;
}



template <typename SCALAR_TYPE>
bool operator != (const class vector_class<SCALAR_TYPE> &V , const class vector_class<SCALAR_TYPE> &W)
{
  return (!(V == W));
}



// Maximal value of array
// ----------------------
template <typename SCALAR_TYPE>
SCALAR_TYPE vector_class<SCALAR_TYPE>::max () const
{
  if (dimension == 0) error_message_print_abort ("max function not defined for vector_classs without elements");

  SCALAR_TYPE List_max = List[0];

  for (unsigned int i = 1 ; i < dimension ; i++)
    {
      const SCALAR_TYPE List_i = List[i];

      if (List_max < List_i) List_max = List_i;
    }

  return List_max;
}

// Minimal value of vector_class
// ----------------------------
template <typename SCALAR_TYPE>
SCALAR_TYPE vector_class<SCALAR_TYPE>::min () const
{
  if (dimension == 0) error_message_print_abort ("min function not defined for vector_classs without elements");

  SCALAR_TYPE List_min = List[0];

  for (unsigned int i = 1 ; i < dimension ; i++)
    {
      const SCALAR_TYPE List_i = List[i];

      if (List_min > List_i) List_min = List_i;
    }

  return List_min;
}



// Memory occupied by a vector in Mb
// ---------------------------------
template <typename SCALAR_TYPE>
double used_memory_calc (const class vector_class<SCALAR_TYPE> &T)
{
  return ((sizeof (class vector_class<SCALAR_TYPE>) + T.get_dimension ()*sizeof (SCALAR_TYPE))/1000000.0);
}

#endif



